Numerical blowup in two-dimensional Boussinesq equations 
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In this paper, we perform a three-stage numerical relay to investigate the finite time singularity in 
the two-dimensional Boussinesq approximation equations. The initial asymmetric condition is the 
middle-stage output of a 2048 run, the highest resolution in our study is 40960 2 , and some signals 
of numerical blowup are observed. 

PACS numbers: 47.20.Cq, 47.27.Te, 47.27.Eq, 47.27.Jv 



It is still an open question whether smooth initial con- 
ditions in three-dimensional (3D) Euler equations can 
develop singularities in a finite time. The current nu- 
merical methods can not really express an infinite num- 
ber in a dynamic simulation, so besides the right initial 
conditions, extremely high resolutions are necessary to 
catch the signal of blowup. The obvious ways to increase 
the resolution are: 1) Using the largest computer avail- 
able to perform simulations; 2) Adoptingsome kinds of 
symmetries in the 3D Euler equations [1, 0, H. The 
two-dimensional (2D) Boussinesq approximation equa- 
tions correspond to the 3D axisymmetric Euler equa- 
tions, and need much less computer capacity than those 
3D symmetric models. In the meanwhile, it can reveal 
much more physics than the one-dimensional symmetric 
model 0, S| • The blowup signal of 2D Boussinesq equa- 
tions is derived in @, 0] : if the maximum absolute values 
of the vorticity and temperature gradient behave like 



(t c -t)- a k (t c - ty 



(1) 



with a > 1 and [3 > 2, a finite time singularity will be 
developed. 

The equations under consideration are the follow- 
ing @: 



e t + u • V0 = o, 

u>t + u • Vw = —t 

A-0 = -u, 



(2) 
(3) 
(4) 



where 9 is the temperature, u = (u,v) the velocity, u) — 
(0, 0, u>) = V x u vorticity, and tp stream function. 

The method adopted in our numerical simulations is 
the filter pseudo-spectral method with some proper de- 
aliasing technique. The modifying factor for each Fourier 
mode k in the filter is ip(k) = e ~ 37< - 2k / N) for k < N/2, 
where N is the grid number in one direction. The de- 
aliasing scheme to be used is the phase-shift scheme with 
circular truncation [8j . 



The low-storage third-order Runge-Kutta time dis- 
cretization is adopted. Besides its high accuracy, the 
one-step property of this scheme is essential to our three- 
stage numerical relay Q. 

During the process towards the singularity, a 5 or 5- 
like function will be developed in the flow field. It is 
well known that spectral coefficients of the 6 function 
are constant for different modes: 



S(x,y)~C 



D —i(nx-\-my) 



(5) 



where C is a constant. If the resolution A 2 — > 00, we will 
get the exact Fourier representation, which is impossible 
in current computers. The main idea of this research is 
to make A" as large as possible to capture the blowup 
signal. 

For high-resolution simulations, it is very important to 
have an effective initial condition because a poorly cho- 
sen initial condition may not lead to blowup or may lead 
to blowup after an unacceptable long-time computation. 
In the following, we will describe how to obtain the initial 
data for the whole paper. First, we take the initial con- 
dition with unified zero vorticity and a cap-like contour 
of temperature with the following expression: 



0(x,y,0)=50( 



Ax — 37T 



7T 



)6 1 (x,y)e 2 (x,y) [1 - 9i(x, y)] , 

(6) 

(x — 7r) 2 is posi- 
and zero otherwise; 
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where if S(x,y) 
tive, 61 = exp (l — tt 2 /S(x, y)) 
if s(y) := \y — 2n\ /1.957T is less than 1, 9 2 = 
exp (l — (1 — s(j/) 2 ) -1 ) , and zero otherwise. We com- 
press the intermediate results (at t = 1.2) starting from 
the above initial condition to form a new initial data. 
More precisely, we let u(x, y, 0) = w'{x, 2y—QA-K, 1.2) and 
9(x,y,0) = 6'(x,2y-0Air,1.2), for (x,y) 6 [0,2tt] x [0,tt] 
(where 9' and u' are obtained by solving Eqs. ©-(|1]) 
and© with a 2048 2 grid), and zero otherwise. To elim- 
inate the high order frequency generated from the com- 
pression, we perform a 2048 2 run with the new initial 
data. The intermediate results at t = 0.12 are the REAL 
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vorticity contours at 1=0 



density contours at t=0 




vorticity contours at t=0.75 



density contours at t=0.75 




FIG. 1: Contour plots of temperature and vorticity with the 
resolution of 8192 2 . Here, "x" indicates the location of uj ma x, 
and "y" the location of |V#| marc . 



initial data for the whole paper. The flow field in spectral 
space is stored, with 2048 modes in both directions. 

The system describes a cap- like hot zone of fluid rising 
from the bottom, while the edges of the cap lag behind, 
forming eye- like vortices (Fig. (JTJ). The hot liquid is 
driven by the buoyancy and meanwhile attracted by the 
vortices, which leads to the singularity-forming mecha- 
nism in our simulation. The flow filed is asymmetric to 
x = 7r to avoid other mechanism during the singularity 
forming 0, @] . Three-stage simulations are planned: 

Stagel Full-time simulations are carried out for five res- 
olutions: 1024 2 , 2048 2 , 4096 2 , 6144 2 , and 8192 2 . 
The time steps are 1.0 xlfr 3 , 6.0 x 10~ 4 , 3.0 xl0~ 4 , 
2.0 x 10 -4 , and 1.5 x 10 -4 , respectively, given by 
the CFL condition. For the 1024 2 run, the first 
1024 modes of all 2048 modes in both directions 
are adopted for initial data. 

For grids finer than 2048 2 , the lower modes use the 
2048 2 initial values, and higher modes are all zero 
(Stage2&3 adopt the similar strategy). 

Stage2 The intermediate results at t = 0.75 of the 8192 2 
run are used as the starting point for three resolu- 
tions: 12288 2 , 16384 2 , and 20480 2 . The time steps 
are 1.0 x 10~ 4 , 7.5 x 10~ 5 , and 5.0 x 10~ 5 , respec- 
tively. 

Stage3 The intermediate results at t = 0.84 of the 
20480 2 run are used as the starting point for three 
resolutions: 24576 2 , 32768 2 , and 40960 2 . The time 
steps are 4.0 x 10~ 5 , 3.75 x 10~ 5 , and 2.5 x 10~ 5 , 
respectively. 
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FIG. 2: The evolution of the T2 and maximum 6 errors for all 
resolutions. The errors are denned as (T2(to) — 2Mi))/T2(io) 
and \9 max {to) - 9 max (t)\/\8 max (t )\. Note that t = for 
Stagel, t = 0.75 for Stage2, and to = 0.84 for Stage3. 



In time-evolution plots of uj max and \V0\ max (Figs. 
® ©), the overlap between 6144 2 and 8192 2 curves lasts 
after t = 0.8, which makes us believe it is safe to use the 
t = 0.75 results of 8192 2 as a starting point for Stage2. 
The starting point for Stage3 is decided similarly. All 
simulations are carried out after the possible blowup time 
with our MPI C++ solver (lb] . 

To demonstrate the accuracy of our numerical results, 
we will check two values which should be time indepen- 
dent during the simulations due to the divergence-free 
constraint, the doubly-periodic condition and the invis- 
cid transport equation (Eq. ©): 



(x,y,t)dxdy, 



Fig. [2] shows that these global average quantities are 
well conserved for all simulations, and errors are better 
controlled for finer grids. For the largest resolution used 
in this paper (40960 2 ), the T 2 error is below 1.0 x 10~ 7 , 
and the 8 max error is below 3.0 x 10 -8 until the supposed 
blowup time is closing. The 1024 2 run has a very poor 
performance according to Fig. O We show it here mainly 
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FIG. 4: Time evolutions of maximum u. 




FIG. 5: Time evolutions of maximum \V6\. 



because it is the largest resolution currently adopted in 
those full 3D Eulcr investigations on this issue. 

If a singularity is about to form at t — T c on the point 
{xcUc), the locations of uj max and |V6*| mQ:! ; should ap- 
proach (x c ,y c ) when t — > T c , or, the distance of the lo- 
cations of these two peak values should go to zero as 
t — > T c . According to Fig. [3J the distances between 
these two peak values in 2048 2 , 4096 2 , 6144 2 , and 8192 2 
runs experience similar drop-down and increasing pro- 
cess until t « 0.846 when they suddenly drop down to 
0.1. The simulations in Stage2&3 reveal some more de- 
tails: it seems that there is another drop-down to almost 
zero at t « 0.86, it seems that there might be a singular- 
ity forming around that time, and we will focus on this 
point in the following. 

Although we eliminate the high order mode distur- 



bance by setting all of them to zero in this study, time 
evolutions of maximum u) in Stagel are still quite similar 
to previous investigations @, Q (Figs. ED. During the 
process of increasing resolutions, the uj max values at the 
late time are always getting higher. The higher resolu- 
tions we adopt, the larger parts of the u> max evolution 
curves of two neighboring resolutions will overlap. For 
example, the overlap between 1024 2 and 2048 2 lasts un- 
til t = 0.18, while that of 2048 2 and 4096 2 lasts until 
t = 0.39, and that of 4096 2 and 8192 2 lasts until t = 0.5. 
The u m ax time evolutions in Stage2&3 are distinguished 
from Stagel. Although higher resolutions still lead to 
higher u> max , it only happens after t 0.86. And just 
before t « 0.86, the values of uj max converge for differ- 
ent resolutions in Stage2&3. So, unlike Stagel, for finer 
grids, the overlap parts between two neighboring curves 
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vorticity contours at t=0.84 density contours at t=0.84 




4.6 4.8 5 5.2 5.4 5.6 4.6 4.8 5 5.2 5.4 5.6 

vorticity contours at t=0.8588 density contours at t=0.8588 
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vorticity contours at t=0.86 density contours at t=0.86 
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FIG. 6: Contour plots of temperature and vorticity near the 
blowup time with the resolution of 40960 2 . Here, "x" indi- 
cates the location of uj ma x, and "y" the location of |V0| ma x. 

will not become longer in Stage2&3. 

From the discussion in the above paragraph, we realize 
that although our numerical scheme is globally adopted 
in Computational Fluid Dynamics (CFD) with solid the- 
oretical proof [9(, it can not generate more accurate re- 
sults for the current problem when the resolution is larger 
than 12288 2 . A common knowledge in CFD field is: 
If the flow field is smooth, a proper numerical 
scheme with higher resolution will lead to more 



accurate results. Its inverse negative proposition is 
also the truth: If a proper numerical scheme with 
higher resolution will not lead to more accurate 
results, the flow field is not smooth. It is clear that 
there is a singularity at t » 0.86 in our problem. 

In the simulation, the values of |V0| are not filtered 
directly. Time evolutions of |V(9| maa; are, by all means, 
the traditional ways: for larger resolutions, the |V0| maa; 
values at the late time are always higher; the higher reso- 
lution we adopt, the larger overlap between the \\70\ m ax 
evolution curves of neighboring resolutions is. Using the 
t e [0.84, 0.859] values of the 40960 2 run, we estimate 
that T c = 0.86, a is slightly larger than 1, and /3 = 2.891 
(Eq. [I]). 

The physical process around t = 0.86 is shown in Fig. 
[6] There is a secondary vortex forming near [5.3, 1.5] 
since t — 0.84, and it eventually becomes the strongest 
vortex in the whole domain at t — 0.86, and absorb both 
maximum lo and |V0| within its region. 

To sum up, there are several points that make us be- 
lieve that there is a singularity at t rj 0.86: 

• The distance between u) max and \ V9\ max suddenly 
drop down to almost zero around t rj 0.86; 

• For enough fine resolutions, the co rnax (t) curves con- 
verge before t w 0.86, and they suddenly diverge 
after t w 0.86; 

• p > 2. 

This work is supported by National Natural Science 
Foundation of China (G10502054) and the Knowledge 
Innovation Program of the Chinese Academy of Sciences 
(Grant No.KJCX2-YW-L08). Simulations of Stagel were 
finished on local Lenovo Deepcom 1800 supercomputer, 
and those of Stage2&3 were carried out on Dawning 
5000A (Magic cube) in Shanghai Supercomputer Center. 



[1] G.I. Taylor & A.E. Green, Mechanism of the production 
of small eddies from large ones, Proc. Roy. Soc. A 151, 
421 (1935). 

[2] M.Brachet, D.I. Meiron, S.A. Orszag, B.G. Nickel, R.H. 
Moft & U. Frisch, Small-scale structure of the Taylor- 
Green vortex, J. Fluid Mech. 130, 411 (1983). 

[3] S. Kida, Three-dimensional periodic flows with high- 
symmetry, J. Phys. Soc. Japan 54, 2132 (1985). 

[4] Z.Yin & T.Tang, Numerical investigations on the fi- 
nite time singularity in two-dimensional Boussinesq equa- 
tions, |ArXiv:physics/0610053"1 (2006). 

[5] Z.Yin & T.Tang, Resolving Small-scale Struc- 
tures in Two-dimensional Boussinesq Convec- 
tion by Spectral Methods with High Resolutions, 
|ArXiv:physics/0509170| vl (2005). 



[6] A.J. Majda and A.L. Bertozzi, Vorticity and incompress- 
ible flow (Cambridge, 2002). 

[7] W. E & C. Shu, Small-scale structures in Boussinesq con- 
vection, Phys. Fluids 6, 49 (1994). 

[8] G.S. Patterson & S.A. Orszag, Spectral calculations of 
isotropic turbulence: Efficient removal of aliasing inter- 
action, Phys. Fluids 14, 2538 (1971). 

[9] C. Canuto, MY. Hussaini, A. Quarteroni, and TA. 
Zang, Spectral methods in fluid dynamics, Springer, New 
York, Berlin, 1987. 
[10] Z. Yin, L. Yuan & T. Tang, A new parallel strategy 
for two-dimensional incompressible flow simulations us- 
ing pseudo-spectral methods, J. Comput. Phys. 210325 
(2005). 



